#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 
#include "system.h"

void Mdloop(void)
{ 
  int i,j;
  double Av[5],Av2[5],tmp,Tempz,sampleNumber;
  FILE *Fileenergy,*FilePtr;
  

  FilePtr=fopen("movie.pdb","w");
  Fileenergy=fopen("energy.dat","w");
  fprintf(Fileenergy,"#  i    Total    Potential    Kinetic \n");
 
  for(i=0;i<5;i++) { 
    Av[i]=0.0;
    Av2[i]=0.0;
  }

  // initialize radial distribution function and mean-square displacement
  SampleRDF(INITIALIZE);
  SampleDiff(INITIALIZE);

  // loop over all cycles
  for(i=0;i<NumberOfSteps;i++) {
	  
    //Obtain the potential energy and the ideal gas part of the pressure
    Force();
    
   //Calculate the Kinetic Energy
    UKinetic = 0;
    for (j=0; j<NumberOfParticles; j++) 
	    UKinetic+=0.5*(SQR(Velocities[j].x)+SQR(Velocities[j].y)+SQR(Velocities[j].z));
    
    if(i<NumberOfInitializationSteps) {
	    tmp=sqrt(Temperature*3.0*NumberOfParticles/(2*UKinetic));
            for(j=0;j<NumberOfParticles;j++)  {
		    Velocities[j].x*=tmp;
		    Velocities[j].y*=tmp;
		    Velocities[j].z*=tmp;
             } 
             UKinetic *= tmp*tmp;
    } else {
            // add the kinetic part of the pressure
            Pressure +=2.0/3*UKinetic/CUBE(Box);
            UTotal = UKinetic+UPotential;
            Tempz = 2.0*UKinetic/(3.0*NumberOfParticles);
            UKinetic /= NumberOfParticles;
            UPotential /= NumberOfParticles;
            UTotal /= NumberOfParticles;
    
            Av[0]+=Tempz;
            Av[1]+=Pressure;
            Av[2]+=UKinetic;
            Av[3]+=UPotential;
            Av[4]+=UTotal;
            Av2[0]+=Tempz*Tempz;
            Av2[1]+=Pressure*Pressure;
            Av2[2]+=UKinetic*UKinetic;
            Av2[3]+=UPotential*UPotential;
            Av2[4]+=UTotal*UTotal;
	    fprintf(Fileenergy,"%4d  %8g   %8g   %8g\n",i-NumberOfInitializationSteps,UTotal,UPotential,UKinetic);
      
           if((i%500)==0)
		    printf("Step: %d UTotal/N: %lf UKinetic/N: %lf UPotential/N: %lf Temperature: %lf Pressure: %lf\n",
                           i,UTotal,UKinetic,UPotential,Tempz,Pressure);
			   
			   
           if(i%50==0) WritePdb(FilePtr);
     
           // sample radial distribution function and mean-square displacement
           SampleRDF(SAMPLE);
           SampleDiff(SAMPLE);
    }
    
     // integrate the equations of motion, get the positios and velocities of next step
     Integrate();
    
  }
 
  // print averages to screen
  sampleNumber=NumberOfSteps-NumberOfInitializationSteps;
  for(i=0;i<5;i++) {
	    Av[i] /= sampleNumber;
    	    Av2[i] /= sampleNumber;
   }
   printf("\n");
  //print the average and standard error
   /*printf("Temperature  : %lf (%lf)\n",Av[0],sqrt((Av2[0]-Av[0]*Av[0])/(sampleNumber-1)));
   printf("Pressure     : %lf (%lf)\n",Av[1],sqrt((Av2[1]-Av[1]*Av[1])/(sampleNumber-1)));
   printf("UKinetic/N   : %lf (%lf)\n",Av[2],sqrt((Av2[2]-Av[2]*Av[2])/(sampleNumber-1)));
   printf("UPotential/N : %lf (%lf)\n",Av[3],sqrt((Av2[3]-Av[3]*Av[3])/(sampleNumber-1)));
   printf("UTotal/N     : %lf (%lf)\n",Av[4],sqrt((Av2[4]-Av[4]*Av[4])/(sampleNumber-1)));*/
   //print the average and standard deviation
   printf("Temperature  : %lf (%lf)\n",Av[0],sqrt(Av2[0]-Av[0]*Av[0]));
   printf("Pressure     : %lf (%lf)\n",Av[1],sqrt(Av2[1]-Av[1]*Av[1]));
   printf("UKinetic/N   : %lf (%lf)\n",Av[2],sqrt(Av2[2]-Av[2]*Av[2]));
   printf("UPotential/N : %lf (%lf)\n",Av[3],sqrt(Av2[3]-Av[3]*Av[3]));
   printf("UTotal/N     : %lf (%lf)\n",Av[4],sqrt(Av2[4]-Av[4]*Av[4]));
  

  SampleRDF(WRITE_RESULTS);
  SampleDiff(WRITE_RESULTS);

  fclose(Fileenergy);
  fclose(FilePtr);
}
